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ABSTRACT 

We performed a numerical experiment designed for core formation in a self-gravitating, 
magnetically supercritical, supersonically turbulent, isothermal cloud. A density prob- 
ability distribution function (PDF) averaged over a converged turbulent state before 
turning self-gravity on is well-fitted with a lognormal distribution. However, after 
turning self-gravity on, the volume fractions of density PDFs at a high density tail, 
compared with the lognormal distribution, increase as time goes on. In order to see 
the effect of self-gravity on core formation rates, we compared the core formation rate 
per free-fall time (CFRg ) from the theory based on the lognormal distribution and the 
one from our numerical experiment. For our fiducial value of a critical density, 100, 
normalised with an initial value, the latter CFRff is about 30 times larger the former 
one. Therefore, self-gravity plays an important role in significantly increasing CFRg. 
This result implies that core (star) formation rates or core (stellar) mass functions pre- 
dicted from theories based on the lognormal density PDF need some modifications. 
Our result of the increased volume fraction of density PDFs after turning self-gravity 
on is consistent with power-law like tails commonly observed at higher ends of visual 
extinction PDFs of active star-forming clouds. 

Key words: ISM: clouds - methods: numerical - MHD - stars: formation - turbu- 
lence. 



1 INTRODUCTION 

Probability distribution functions (PDFs) for volume and 
column density fields calculated from turbulence simulations 
with the isothermal equation of state ha ve been successfully 
fitted with a lognormal d ist ribution (IVazquez-Semadenil 
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1998 



Padoan et al.l Il997 : iPassot fc Vazquez-Semaden: 
Nordlund fc Padoanl Il999 |: lOstriker et all ll99E 
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Lemaster fc Stonel 120081 ; iFederrath et all l2008aT ). In addi- 
tion to these numerical works, it has been recently reported 
that the Ha emi ssion measure for the warm ionized medium 
(|Hill et al.l l2008) and the densities of the diffuse ionized gas 
and diffuse atomic gas l|Berkhuiisen fc Fletcher! l200Sh fol- 
low the lognormal distribution. The lognormal density PDF 
has also become one of important ingredients of star for- 
m ation theories based on turb ulenc e. The analytical models 
of iPadoan fc Nordlundl |2002l ) and iHennebelle fc Chabrierl 
(2008) predicted stellar initial or core mass functions, whose 
mass distributions at lower masses are largely determined 
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by the lognormal density PDF. lElmegreenl (|2002l . l2008h 
used the lognormal density PDF in order to calculate the 
cumulative mass fraction above a critical number density 
(> 10 s c m" 3 ), which is directly relate d to the star formation 
efficiency. iKrumholz fc McKeel (|2005l ) also measured dimen- 
sionless core (star) formation rate per free-fall time (CFRh) 
based on the lognormal density PDF. Ins tead of pick i ng up 
a spe cific value of the critical density as in Elmegreen (2002, 
2008), they determined it by equating the local Jeans and 
sonic lengths. 

The numerical simulations mentioned in the first para- 
graph didn't include self-gravity of gas. It is, however, eas- 
ily expected that the inclusion of self-gravity in a turbulent 
simulation results in a density PDF with an extended tail 
at higher densities due to the very nature of self-g ravity, 
which has been, in fact, reported in a few literatures. iKlesser] 
(2000) showed the extended tails of density PDFs from 
his dec ay and driven hydrody namic simulations with self- 
gravity. iDib fc Burkertl (|2005l ) reported a similar extended 
tail of a density PDF from one of thei r numerical models 
with self-gravity. IFederrath et alj (|2008al ) showed that power 
laws develop at the high density tail of the density PDFs 
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Figure 1. Volume renderings of two density fields at t/tgf = 0.0 (left) and t/tg = 1.2 (right), respectively, where iff is the free-fall time 
with an initial density. The evolutionary time is set to zero when self-gravity is on. Colours in a bar are mapped onto density values 
normalised with an initial value, po- 



of tracer particles in self-gravitating, supersoni c turbu lence 
calculations. Finally, IVazquez-Semadeni et all (|2008T ) per- 
formed self-gravitating, isothermal, hydrodynamic simula- 
tions and showed progressive increase of the volume fraction 
of density PDFs at high densities as a function of time. 

The motivation of this paper is to show quantitatively 
how much CFRff in a turbulent cloud can be increased by 
self-gravity with respective to that measured based on a log- 
normal density PDF. As we mentioned above, density PDFs 
from isothermal simulations without self-gravity result in 
log-normal functions. Dense filaments formed in those simu- 
lations are transient, which reside in the high end of the func- 
tions and are the main contributor to core or star formation 
rates. It is self- gravity that enables some of them to develop 
into collapsing cores, and helps the cores to accrete nearby 
gas at later evolutionary stages. With this motivation in our 
mind, we perform a numerical simulation for core formation 
in a self-gravitating, magnetically supercritical, supersoni- 
cally turbulent, isothermal cloud. We show that the volume 
fractions of density PDFs at higher densities are increasing 
as time goes on. Even though the increments of the volume 
fractions at high densities over a lognormal distribution are 
small, the increments of mass fraction are quite large due to 
the very high densities. This, in fact, leads to more signif- 
icant increment of CFRff than the one measured based on 
the lognormal density PDF. 



2 NUMERICAL METHOD 

In order to see the formation and evolution of cores in a tur- 
bulent cloud, we numerically integrate the MHD equations 
with an isothermal equation of state and the Poisson's equa- 
tion, using a MHD code based on a total variation dimin- 
ishing scheme l|Kim et al.lll999l ) and a gravity solver based 
on the fast Fourier transform. As an initial condition, a uni- 
formly magnetized medium is assumed in a rectangular box. 
The periodic boundary condition is imposed in each direc- 
tion of the box. The combined isothermal MHD and Pois- 
son's equations are scale-free, and can be written in a di- 
mensionless form with two parameters. One is the plasma 
beta, /3, the ratio of gas to magnetic pressures, and the other 



is the Jeans number, J, the ratio of the length of one side 
of the computational box, L, to the initial Jeans length, 
Ajo = i^oi/Gpo) 1 ^ 2 , where a g is an isothermal sound speed, 
G is the gravitational constant, and po is an initial density. 

We basica l ly follo w th e recipes for turb ulence generation 
in lStone et all (l99Sl ) and lMac Lo"wl (|l999l ). Velocity fluctu- 
ations are generated in a Fourier space with th e same func- 
tional form of the velocity power spectrum in IStone et al.l 
(1998). We, however, take the peak wavenumber, 2(2iv/L), 
which is smaller than their choice. We then transform the 
fluctuations into a real space, and adjust their amplitudes 
in a way that an input kinetic energy rate is a constant. 
The level of turbulence driven by this method can be pa- 
rameterized by a rms (root-mean-square) sonic Mach num- 
ber, M a . We start a driven turbulence simulation without 
turning self-gravity on, and wait until a converged turbu- 
lent flow is developed. Then, we turn self-gravity on and 
set the time to zero. The forcing for turbulence generation 
is still active afterwards. Numerical simulations based on 
the above-mentioned method with different parameters have 
been done. In this Letter, we show results from one 512 3 
simulation with parameters, /3 — 0.1, J = 4, and M s = 10. 



3 RESULTS 

To visualise the evolution of density fields as a function 
of time, the series of volume rendering images of three- 
dimensional density fields are made. Two of them are shown 
in Figure 1. The left image shows a density field at time 
just before turning self-gravity on, t/tg = 0, where tg = 
[37r/(32Gpo)] 1 ^ 2 , the free-fall time with an initial or a mean 
density inside the computational box. The colour bar only 
covers a normalised density range from 1 to 2000, which 
enables us to see the distribution of high density gas more 
clearly. There are many dense filaments in the left image. 
These are formed by the interaction of large-scale super- 
sonic flows with a sonic Mach number around 10. Isothermal 
shocks along a field line easily increase the normalised den- 
sities at post shock regions up to around 100, which is the 
square of the Mach number. After turning self-gravity on, 
some of the filaments increase their central densities either 
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Figure 2. Evolution of the volume fractions of density PDFs. 
The horizontal and vertical axes are logarithmic values based on 
e and 10, respectively. A black line in each panel is a density PDF 
in a computational box at a specified time in units of the initial 
free-fall time. A red line is a lognormal distribution, equation 1(1)1 . 
fitted with a density PDF averaged over the saturated stage of 
turbulence before turning self-gravity on. The fitted parameters 
are /x = —0.82 and a = 1.28. The same red line is plotted in each 
panel. 



by merging with each other or accreting nearby gas. Several 
condensations are seen in red or yellow colours in the right 
image at t/tg = 1.2, whose normalised central densities are 
larger than 300. It is these condensations that are eventu- 
ally evolved into collapsing cores. In fact, they are the main 
contributors to form bumps raised on top of a lognormal 
density PDF at high densities (see Figure 2). 

We plot two kinds of density PDFs in black and red 
lines in each panel of Figure 2. The black lines are calculated 
from density fields from our numerical simulation at different 
times. A same red line in each panel is a fit of an averaged 
density PDF over a saturated turbulent state before turning 
self-gravity on with a lognormal distribution, 



Pln(s) = 



V2tt<t 2 



exp 



2a 2 



(1) 



where s = ln(p/po), M is t ne mean, and a is the stan- 
dard deviation. Since mass has not added into or sub- 
tracted from the computational box during our simula- 
tion, total mass inside the box is conserved. The conser- 
vation constraint at the initial state and a state where the 
density PDF is described with the lognormal distribution, 




Figure 3. Cumulative mass fractions as a function of the nor- 
malised density. Solid lines with different colours are calculated 
by integrating equation J2j numerically with the density PDFs 
at different times shown in Figure 2, respectively. A dash line is 
for the lognormal distribution, which is from equation with 
a = 1.28. 



po 



= ppLN(s)ds, provides us a relation, 



A* = 



-/2 



|Passot fc Vazquez-Semadenil [l998h . This reduces the two 
parameters in equation |l} into one. The average for the 
red line is taken over from t/tg = —1.3 to t/tg = 0.0. The 
number of density fields used in the average is 50. The mean 
value of the averaged density PDF, p(s), is -0.82 calculated 
from "EiSip(si), where Si is a discrete value of s, and p(si) 
is a volume fraction at Si. A lognormal distribution with 
a — 1.28 calculated from the relation a = (— 2/i) 1 ^ 2 is plot- 
ted in each panel as a red line. It serves as a fiducial line to 
see how much a black line at each panel deviates from the 
log-normal distribution. Later, we will also use the a value 
to estimate CFRjf based on the lognormal distribution. 

A black line in the upper left panel at t/tg — shows a 
density PDF at the turn-on time of self-gravity. It is quite 
well matched with the red line, except at low and high 
density ends. The differences may be due to finite numer- 
ical r esolution and the intermittency effect of t he turbu- 
lence dKritsuk et al ] |2007l ; lFederrath et al. I l2010al ). We then 
deliberately choose several times based on shapes of den- 
sity PDFs. The shape of density PDFs, as time goes on, 
is hardly changed from the log-normal distribution up to 
around t/tg — 0.91. It takes for self-gravity to exercise its 
control over supersonic turbulent flows and then bring sig- 
nificant change in density PDFs. In fact, it is at t /tg = 1.05 
that the black line shows quite an excess of volume fractions 
at a high density tail, which we call a bump on top of the log- 
normal distribution. The bump height is increasing rapidly 
up to t/tg = 1.44, and then almost saturated afterwards. 
Even though the excess volume fractions over the lognor- 
mal distribution is small, the excess mass fractions should 
be quite significant due to high densities at the bump. It 
is this excess mass fractions that contribute significantly to 
core formation rates calculated in the followings. 

We calculate the cumulative mass fractions of the den- 
sity fields. Since the volume fraction of a density PDF at 
each time is known, the cumulative mass fraction, f(s), of 
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Figure 4. One-dimensional velocity dispersion as a function of 
box size. Normalization units for the dispersion and the size are 
the isothermal speed a s and the one-dimensional size L of the 
computational box. respectively. For a given time, measured mean 
values of velocity dispersions at different box sizes are connected, 
and the standard deviations are plotted with error bars. The nor- 
malised times t/tfj and sonic lengths X a /L are given in parenthe- 
ses. 



gas whose density is larger than s is easily calculated using 
the following equation, 



/(*)= / -P(s)ds 
Po 



e s p(s)ds. 



(2) 



If p(s) is replaced with equation |T| with the relation [i — 
— (j 2 /2, the integration results in 



/ln(s) = - 



1 + erf 



cr 2 - 2s 

2\/2cr 



(3) 



which is the same as the expression in equation (20) for 
the core formation r ate, i .e., ignoring feedback effects, in 
iKrumholz fc McKed l|2005h . We numerically integrate equa- 
tion ((2| with p(s)'s shown in Figure 2, and then plot cu- 
mulative mass fraction curves as a function of log(p/po) in 
Figure 3. A dash line is calculated from equation @ with 
g = 1.28. A red line at t/ts = 0.0 always lies below the 
dashed line. This is due to the fact that the density PDF at 
t/tfj = 0.0 drawn with a black line shown in Figure 2 has 
smaller volume fractions at high densities than the red line. 
As time goes on, the cumulative mass fraction increases sig- 
nificantly at high densities. As it is expected from Figure 2, 
the increasing rate is high during the rapid development of 
the bump from t/ts = 1-05 to 1.44. Furthermore, the mass 
fraction curves at later times are significant higher than the 
one calculated from the lognormal distribution, which is to- 
tally due to self-gravity. 

The cumulative mass fraction curves shown in Fig- 
ure 3 may provide CFRff, if a critical density for gas col- 
lapse is properly defined. One idea that has been used is 
to compare a local Jeans length , Aj = (Tia^/Gp) 1 ^ 2 , and 



a sonic length, A., JPadoan 19951 ; IVazquez-Semadeni et al.l 



2003; IKrumholz fc McKeel 12005^ In fact, the condition of 
Aj = A s gives a critical density p c normalised with the ini- 
tial density, 



El - ( ^ J0 \ 
Po I A, / 



(4) 



where Ajo is the initial Jeans length. The physical back- 
ground of this condition is that a transonic turbulent ve- 
locity dispersion over a volume defined by A s is barley able 
to support gravitational collapse of the gas inside the vol- 
ume. Since Ajo = L/4, which is given as an initial condi- 
tion, we need X 3 in te rm of L. For the calculati o n of A s , we 
follow the method in IVazquez -Semade ni et al. | i|2003h . We 
take 100 random positions in the computational box, put a 
same-sized box less than L centered at each position, and 
calculate one-dimensional velocity dispersions of turbulent 
velocity fields in the 100 boxes. Figure 4 shows the velocity 
dispersions normalised with the isothermal sound speed as 
a function of a normalised box size. The dispersions at eight 
different times are plotted with coloured solid lines. Error 
bars at the measured box size are included, even though 
they are hardly distinguished from each other. A horizontal 
dash-dot line with a unit normalised velocity dispersion is 
drawn. In fact, sonic lengths are determined by the meeting 
points of the dash-dot line and solid lines. The measured 
sonic lengths at different times are included in the panel, 
which are in the range, 0.033L 5 A s ^ 0.050L. If this range 
of A s values and Ajo = L/4 are plugged in equation @, 
then 25 ^ pc/ Po ~ 57. Considering the fact that a density 
jump brought by a Mach 10 isothermal shock along a mag- 
netic field line is around 100, the range of the normalised 
critical density can be easily attained by shocks only. So it 
is not guaranteed for the gas to collapse, whose normalised 
density is the upper bound of the critical density range, 57. 
Furthermore, the magnetic field included in our simulation 
plays a certain role in supporting cores. So the argument 
based on the local Jeans and sonic lengths may not give a 
correct critical density for measuring CFRff in our numerical 
simulation. 

In order to truly measure CFRff in our simulation, we 
first calculate core formation efficiency, which is the ratio of 
the total mass of cores defined by a critical density to the 
total mass in the computational box. Figure 5 shows the 
evolution of core formation efficiencies with several different 
critical densities. Because of the reason given in the pre- 
vious paragraph, we choose a rather wide range of critical 
densities normalised with the initial density, from 30 to 500. 
The core formation efficiencies shown in red and brown lines 
show more or less constant levels up to around t — 0.8tg. 
There is very tiny fraction of the total mass during the inter- 
val, whose normalised density is larger than 100. This is due 
to the fact that it takes for self-gravity to exercise its con- 
trol over turbulent flows, as we have seen in the PDF plots 
(Figure 2). After that point, the efficiencies increase rapidly 
up to around t = 1.5tg and then more slowly later on. The 
increase of the efficiencies is mostly due to the accretion of 
nearby gas onto several cores. 

For a given critical density, CFRff can be calculated 
using equation (|3}. Here the free-fall time is again mea- 
sured based on a mean density of a molecular cloud. For 
the seven normalised critical density values, 30, 50, 100, 
200, 300, 400, and 500 shown in Figure 5, the CFR ff 's are 
2.2 x 10~ 2 , 7.8 x 10" 3 , 1.5 x 10" 3 , 2.3 x 10~ 4 , 6.8 x 10" 5 , 
2.7 x 10 _J , 1.2 x 10 -5 , respectively. In order to measure the 
CFRff's from our numerical simulation, we make a least- 
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Figure 5. Evolution of core formation efficiencies with different 
critical densities. The core formation efficiency is defined as a 
fraction of total mass in the computational box, whose density 
is larger than a critical density. The time is measured in units of 
free-fall time with an initial density. The critical density, p c , is 
normalised with the initial density, pg. 



square fitting of each curve from t = to t = 3tft in Figure 
5 with a straight line. In fact, the slopes of those lines give 
us CFRff's. They are 0.048, 0.045, 0.040, 0.036, 0.033, 0.031, 
0.029 for normalised critical density values, 30, 50, 100, 200, 
300, 400, and 500, respectively. Comparison of theoretical 
and numerical estimates shows that the CFRfi's from the 
numerical experiment are 2.2, 5.8, 27, 160, 490, 1100, and 
2400 times larger than those from the lognormal distribu- 
tion for the critical density values 30, 50, 100, 200, 300, 400, 
and 500, respectively. The difference becomes larger as the 
critical density increases. If we pick 100 up as our fiducial, 
normalised critical density, the CFRff based on the theory 
is likely to be underestimated by about a factor 30. 



4 DISCUSSIONS AND CONCLUSIONS 

Recently iKainulainen et alj <|200 9) catalogued column den- 
sity PDFs of 23 molecular cloud complexes from the 2MASS 
archive. They classified them into two groups based on star 
formation activity and compared their column density PDFs 
with each other. The column density PDFs of star-forming 
clouds always have extended tails, whereas the PDFs of 
clouds without active star formation follow lognormal dis- 
tributions or a bit excess at high column densities. Further- 
more, the cumulative fractions of column density PDFs with 
star-forming clouds are significantly larger than those with- 
out active star formation. These observational results clearly 
show that self-gravity plays a role in forming the extended 
tails of the PDFs, which is consistent with our work. 

We remind that the extended tails of density PDFs at 
high densities from isothermal simulations with self-gravity 
have been shown in a few previous literatures. The new 
finding in this Letter is that the extended tails can en- 
hance CFRff quite significantly. Cores in a turbulent cloud 
cannot form without self-gravity. Core (star) formation 
rates or core (stellar) mass functions should be measured 



in the context of a self-gravitating cloud. Therefore, it is 
likely that core formation rates previously measured based 
on the lognorma l density PDF jKrumholz fc McKeel 120051 : 
lElmegreenll2008l ) are underestimated. Likewise, the core or 
stellar i nitial mass functions based on the lognormal distri- 
butio n (jPadoan fc Nordlu nd 2002:; iHennebelle fc Chabrierl 
2008) need to be modified. 

There are uncertainties in our results due to numerical 
resolution, measuring CFRff based on only a density thresh- 
old, violation of the numerical Jeans condition, and stellar 
feedback. Because of the limited space, we briefly discuss on 
the last three. Firstly, the density threshold alone may not 
fully capture collapsing gas. Mor e elaborated coll a pse ind i- 
cators are needed (for example, iFederrath et al.1 (|2010bh ). 
Secondly, Figure 2, for example, shows that density values 
at a high density tail especially at later stages of our sim- 
ulation go above the ma ximum normalised density valu e, 
1024, (see, Equation (9) in lVazquez-Semadeni et aL ([2005)) 
constrained by the nume rical Jeans condition for preventing 
artificial fragmentation |Truelove et al.|[l997T ). In this Let- 
ter we are interested in not the fragmentation of cores but 
the total amount of mass of cores defined by a critical den- 
sity. So the violation of the condition will not change our 
main result but add uncertainty in the measured CFRff. 
Thirdly, we didn't include feedback processes from the stars 
that might form in our simulation. Without the stellar feed- 
back core formation efficiency eventually approaches one. 
However, at least, before the formation of a first star in the 
simulation, our measurement of the core formation efficiency 
is quite right. In order to properly measure the core forma- 
tion efficiency, especially, at the later evolutionary state of 
a molecular cloud, one should include the feedback. 

We performed a magnetically supercritical, supersonic 
turbulence simulation with the isothermal equation of state 
to study the effects of self-gravity on density PDFs and the 
core formation rate. Here are conclusions from the study. 
First, self-gravity helps to form the extended tail of a density 
PDF at high densities, which significantly increases CFRff. 
Second, the normalised critical density for core collapse de- 
termined by the equal condition between the local Jeans and 
sonic lengths is 25 5; pc/ po ~ 57, which is smaller than 100, 
a density jump brought by an isothermal Mach 10 shock in 
our simulation. So the determined critical density may not 
give a correct condition for the core formation. Third, for 
our fiducial normalised critical density, 100, CFRff = 0.045, 
measured from our numerical simulation is about 30 times 
larger than the one, 0.0015, based on the lognormal dis- 
tribution. Therefore, self-gravity plays a significant role in 
enhancing CFRff in a turbulent cloud. 
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